Bosonic enhancement of spontaneous emission near an interface 
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Abstract 

We show how the spontaneous emission rate of an excited two-level atom placed in a trapped Bose-Einstein condensate of 
ground-state atoms is enhanced by bosonic stimulation. This stimulation depends on the overlap of the excited matter-wave 
packet with the macroscopically occupied condensate wave function, and provides a probe of the spatial coherence of the 
Bose gas. The effect can be used to amplify the distance-dependent decay rate of an excited atom near an interface. 



^—^ L Introduction 

^Spontaneous emission from an excited atom can only take 
^>lace when a vacuum mode is available to accommodate the 
K-^emitted photon. This fact is employed in experiments in cav- 
^ity quantum electrodynamics (cQED) to shift the sponta- 
O^heous emission rate in small cavities where the structure of 
^^;he electromagnetic (em) vacuum is modified PJ [2] . In par- 

! tjcular, it is known that the lifetime of an excited atom near 

JtA, plane surface (the simplest 'cavity' system) shows an oscil- 
l atory behavior for atom-surface distances comparable to the 
-Resonant photon wavelength. The effect is not easy to observe 
^vith ultracold atoms since distance control in the sub-micron 
"J^ange is required. A transient signal related to a change in 
□(lecay rate was observed in the Amsterdam group (3], with 
jcold atoms being probed spectroscopically in the vicinity of a 
' Surface using an evanescent light field. Quite analogous to the 
^ffect of a cavity on photon modes, the presence of a Bose- 
r"Einstein condensate (BEC) can alter the decay of bosonic 
i— -^(toms, as the macroscopic population of atomic modes stim- 
ulates the transition into these HI [51 E] . This enhancement is 



^significantly reduced, however, in a uniform system because 
0£pf momentum conservation. For an excited atom initially at 
£Njest, the final state is shifted by the photon recoil momen- 
'xftum, and overlap with the condensate mode occurs only if 
"Repulsive interactions deplete the ground state [7], or alter- 
^Hiatively, for a confined system where the BEC is spread over 
finite width in momentum space. 

In the present paper, we discuss the enhancement of spon- 
. taneous emission in a trapped BEC and show in particular 
. ^ na t the small oscillations in the decay rate near an surface 
K/ban be significantly amplified. The trap confinement and the 
temperature of the BEC are taken into account and provide 
C^nly a moderate reduction compared to the scaling with the 
number of atoms in the BEC. We use a quantum field theory 
of atoms and photons to calculate complex self-energies (T- 
matrix elements) [8] . The spontaneous emission rate of an ex- 
cited initial state is proportional to the imaginary part of the 
T-matrix element, while its real part yields the Lamb-shift 
of the initial state, and, near a surface, the Casimir-Polder 
potential between BEC and surface, as calculated in [HI EU] 
by different methods. We also apply the formalism to a typ- 
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ical spectroscopy experiment where the enhanced decay rate 
appears in the absorption spectrum of a weak, near-resonant 
laser field. 



2 Decay of an excited wavepacket in a spherical 
atom cloud 

2.1 Self energy and transition rate 

Consider a factorized initial state |e, N) with one atom in the 
electronically excited state and N atoms in the (collective) 
ground state. We apply second-order perturbation theory 
in the interaction with the em field to get for this state a 
complex energy shift (self-energy) whose imaginary part gives 
the Bose-enhanced decay rate. (This formalism is based on 
Refs. |lll [T2"] and explained in detail in Ref.[5].) The em self- 
energy is proportional to the S'-matrix element 

(e,iV|5 (2) |iV,e) = -^ r^^ ^-- (1) 
= J d^xi d i x 2 6(i 2 - h) ( T {E a {x 1 )E j3 {x 2 )} ) 

x (e\^t(x 2 )^ e (x 1 )\e)(N\^ g (x 2 )^l(x 1 )\N} , (2) 

where fi a (a — x, y, z) are the matrix elements of the tran- 
sition dipole. The Feynman diagram represented above uses 
bold lines for the many-body system of ground state atoms, 
dashed lines for individual excited atoms and wavy lines for 
the photon propagator. The brackets (...) in the first line of 
eqn. ^ denote an expectation value with respect to an equi- 
librium state of the em-field, and the symbol T{. . . } denotes 
time-ordering. \P e and ^ g are interaction-picture field oper- 
ators for the atomic levels, with a time dependence governed 
by the unperturbed atomic Hamiltonian. 

In eqn. ([2]), the time dependent phase of the ground state 
correlation function (N\^ g (x 2 )' < i> g (xi)\N} is of the order of 
typical single particle energies in the trap and thus much 
smaller than the phase of the term (e\'$l(x 2 )' l & e (xi)|e) , which 
is proportional to exp[iuj eg (t 2 — t±)] where w eg is the Bohr 
transition frequency. We will hence neglect the time depen- 
dence of the ground state correlation function in eqn. ([2]). 
(This is equivalent to summing the decay rate over the final 
ground state modes.) The time integrations in eqn. ([2| can 



then be performed, yielding for the transition rate (h = 1) 
7 e (7V) = 2[1 + n(u} eg )]fjt a fip J d 3 x |$(x)| 2 ImG tt) g(x,x,a; e9 ) 
+ J d 3 x 1 d 3 x 2 $(x 1 )$*(x 2 )(iV|*t(x 1 )* 3 (x 2 )|iV) 



x ImG Q( g(xi,x 2 , uj. 



eg) 



(3) 



where the excited state is given by the normalized wave func- 
tion $(x). We denote n(ui eg ) the average photon number at 
frequency u) eg in thermal equilibrium. For the sake of sim- 
plicity, we restrict ourselves for the rest of the paper to a 
field at zero temperature where the thermal photon number 
n(w eg ) is negligible. We have expressed the photon propaga- 
tor through the retarded response (Green) function 



G a p{xi,x 2 ) 



did 

2^ 



Xl,X 2 , LO 



= i([E a (x 1 ),Ep(x2)])e(t 1 -t 2 ) 



(4) 



This quantity is easily calculated in a general environment, 
for example near a surface, see Ref.[13], and shows oscillating 
behavior as function of the atom-surface separation. 

In eqn. ([3]), the decay rate naturally splits into a term 
7e°^ that remains in the absence of the atom cloud (first line) 
and an additional term 7 ( P EC ( J ?V) which describes the bosonic 
enhancement. The first term has a natural interpretation in 
terms of an average of local decay rates over the position 
distribution of the excited state wavepacket. (For a study of 
the dynamics of the excited state, see Ref.|14j.) The analysis 
of the second term is the main focus of this paper. Note 
that it depends on the two-point correlation function of the 
ground state atoms. The bosonic stimulation is thus a probe 
of the spatial coherence of the Bose gas. 

2.2 Discussion of Bose enhancement in free space 

To illustrate the physics in eqn. we will assume that both 
the ground-state BEC and the excited atom are trapped in 
overlapping isotropic harmonic potentials, far enough away 
from any macroscopic body. We use the free-space expression 
G ®1 of the em Green tensor (|4|) which takes the form 



G^(xi,x 2 ,o; es ) 



d 3 k 
(2tt) 3 



ik.(xi-x 2 ) 



(5) 



As is well known, only photons on the light cone contribute 
to the the imaginary part of this quantity, i.e., |k| = u) eg /c = 
k eg . Another relevant length scale is the oscillator length 
a = {Mlut)^ 1 ^ 2 of the ground-state trapping potential (M is 
the atomic mass and ujt the trap frequency). We assume that 
the excited-state wave packet ^(x; rf) is an isotropic Gaussian 
with a width rjaQ. The first term in eqn. ^ gives 7e°' — 
/i 2 fc 3 ff /(37T£ ), the free space decay rate [T5] . 

Fig. [I] shows how the decay rate r y e (N) varies with the 
width parameter 77, in the presence of a Bose condensate with 
N = 10 000 Rb atoms (three-dimensional, isotropic trap). 
We consider both an ideal gas model (solid red lines) and 
an interacting Bose gas (dash-dotted lines). The upper black 



curve corresponds to an ideal gas at zero temperature, where 
all the ground-state atoms populate the trap ground state 
■00 (x) , a gaussian with width a,Q . In this case the integrations 
in eqn. ([3| can be worked out explicitly. We find analytically 
that the optimum value of 7 e (-/V) is obtained for 

Vo P t = -y= (y / 9+(fc e5 a ) 4 - (fc^ao) 2 ) 1 / 2 , (0 < 77 opt < 1) . 

(6) 

This result as well as the typical behaviour of j e (N) is easy to 
understand by noting that the decay rate for a given photon 
momentum k is proportional to the overlap integral 



/ 



d 3 a;$(x,7 ? )V;(x) e ik 



(7) 



where n = for the BEC ground mode. In the so-called 
Lamb-Dicke limit k eg ao <C 1 (well- localized trap), the expo- 
nential in eqn. ^ can be approximated by unity, and the 
overlap is optimal when the two wavepackets are matched in 
width, 77 = 1. The opposite case looks closer to a homoge- 
neous system and is easier to analyze in Fourier space where 
the photon recoil provides a shift of the momentum distribu- 
tion. This reduces the overlap and can be compensated for 
by making an excited-state wave packet wider in momentum 
space, i.e., rj < 1. At the optimum value, the width is of 
the order of the photon momentum and the shifted excited- 
state wave packet still has some overlap with the sharp zero- 
momentum component of the BEC. 

The temperature dependence in Fig. [T] closely follows the 
occupation of the ground-state (condensate) mode. We 
have used Ref. |T5] where the two-point correlation function 
(JV|*t(xi)* s (x 2 )|iV) for the ideal Bose gas in a 3D trap is 
given in a simple form, involving only a single summation. 
This provides the Bose enhancement of j e (N) in a straight- 
forward manner for atom temperatures Ta below and above 
the critical temperature T c (see caption). We see that for 
any 77, the transition rate drops below the zero-temperature 
value, and for temperatures above T c , it becomes compara- 
ble to 7^ (horizontal dashed line). The temperature depen- 
dence is shown in fig. [2] and compared to the condensate frac- 
tion Nq/N (dashed line). The bosonic enhancement closely 
follows the population of the condensate mode because the 
excited-state wave packet ^(x) has the largest overlap with 
the trap ground state. The thermal occupation of higher ly- 
ing trap states hence diminishes the integrals eqn. Q. We 
note that this behaviour would change qualitatively in lower- 
dimensional systems where the Bose gas occupies excited 
states with a relatively larger weight. 

The more realistic case of an interacting Bose gas is also 
shown in Figjl] We focus here on repulsive interactions (cor- 
responding to positive s-wave scattering length a s ), and re- 
strict ourselves to temperatures far below T c , where it is le- 
gitimate to approximate the field operator by the condensate 
mode only: 



(JV|*J( X1 )* ( Xa )|JV) a 7VV^(x 1 )0 o (x 2 ) 



(8) 



Elementary excitations of the condensate can be included 
within Bogoliubov theory [HE [T71 \TE\ . The condensate wave 
function -ipo (x) is a solution of the Gross-Pitaevskii equation. 
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We have used the approximate variational solution [T!5] 
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Figure 1: Decay rate of an excited wave packet embedded in a Bose 
condensate with N atoms. The rate "f e (N) [eqn. |3]l] is normalized to 
the free-space value 7e°' = /i 2 fc| g /(37r£o) and plotted as a function of 
the width n arj of the excited state, scaled to the oscillator length for 
the BEC trap. The ground-state cloud consists of N = 10 000 atoms, 
with trap frequency uit/2tt = 1kHz (oscillator length ao = 0.34 /mi), 
the resonance wavelength is 2n/k eg = 780 nm as for rubidium. The op- 
timal value of r] [eqn. |6|] is r? opt = 0.44. Upper black curve: ideal Bose 
gas at zero temperature. Solid (red) curves: ideal Bose gas at tem- 
perature T^4 = 0.3, 0.5, 0.6, 0.8, 0.9, 1.2 T c with critical temperature 
T c = ujt (N I C[3]) x / 3 (top to bottom). Dash-dotted curves: interact- 
ing Bose gas with mode function eqn. (9), for varying s-wave scattering 
length a s = 1, 5, 10a 3 (top to bottom), with k eg a B = 0.047 as for 
rubidium. Horizontal dashed line: free-space decay rate 7e . 
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Figure 2: Decay rate and absorption line width vs temperature Ta of 
an (ideal) Bose gas, for an optimized excited state wave packet. The 
rate j e (N) is given in units of the free-space value 7!°^. Trap and 
atom parameters are the same as in fig. [I] Red (upper) curve: decay 
rate eqn. |3| with n = r) opt = 0.44. Blue (lower) curve: line width 
7abs(^0 [eqn. ( |17| ] of the absorption spectrum. Dashed-dotted curves: 
condensate fraction Nq/N, scaled to the T = values. Horizontal 

dashed line: free-space decay rate 7^. 



MX) = Wrf^ 2 (1 - x2 / R2 ) (1+T) ' 2 @ ( R - x ) < ( 9 ) 

that interpolates between a Gaussian and the Thomas-Fermi 
limit as the parameter Noa s /ao is changing from zero to infin- 
ity. The length R and the exponent r are fixed by minimizing 
the Gross-Pitaevskii energy functional, and Cq is a normaliza- 
tion constant. The result for the decay rate 7 e (JV) is shown by 
the dashed-dotted curves in fig. [lj as the interaction param- 
eter Noa s /ao is increased. Relative to the s-wave scattering 
length a s of rubidium, we took a s = 1, 5, 10a s (top to bot- 
tom) which can be achieved using a Feshbach resonance, for 
example. The interacting gas shows a flatter density profile 
in the trap, as is well known; this results in smaller values of 
the overlap integrals eqn. 

To summarize the data of Figjl] we find a relatively strong 
enhancement of the spontaneous decay rate of an excited 
atom embedded in a Bose condensate. This happens despite 
the non-perfect overlap that encodes the constraints of mo- 
mentum conservation and photon recoil. The optimum condi- 
tions correspond to a well-localized excited-state wavepacket 
(on the scale of the transition wavelength) and a strong con- 
densate fraction (T < 0.5 T c ). 

3 Bose enhancement near a surface 

In this section, we calculate the transition rate eqn. ([3]) near 
a surface and demonstrate its enhancement in a Bose conden- 
sate of oblate shape. This scenario can be realized with an 
optical lattice, by retro-reflecting an off-resonant laser beam 
at the surface [20], or in a bichromatic evanescent wave [21] . 
We take the surface in the xy-plane and the trapped atoms 
centered at a distance d from the surface in the positive z- 
direction. Concerning the surface material, we use the ideal- 
ized model of a perfectly reflecting mirror for the sake of 
simplicity; but with the appropriate choice of (frequency- 
dependent) reflection coefficients that appear in the photonic 
Green function G Q( 3 (xi , x 2 , uS) , a wide range of surface mate- 
rials can be treated in the same maimer [2"2"| 113] . 

As we have measurements in mind where the control over 
the distance d is essential, we take an oblate Bose condensate 
and assume for simplicity a single gaussian mode with widths 
a (in the xy-plane) and a /V\ <C d (along the z-axis). The 
depletion of the condensate and its broadening due to repul- 
sive interactions could be incorporated as in Sec[2j For the 
excited state, we adopt again a gaussian wave packet localized 
in the cloud center, with widths 7700 and ri z ao/VX, respec- 
tively. The actual values of the trap frequency are given in 
the caption of fig. [3] As the very narrow confinement in the z- 
direction describes a quasi-2D scenario, the temperature has 
to be lower than Tj 2 ^ = oju (A r /C[2]) 1/ ' 2 to ensure a strong 
condensate occupation. 

Fig. [3] illustrates the decay rate j e (N) obtained from 
eqn. ^ as a function of the distance d. As is well known, the 
rate depends on the orientation of the dipole moment (par- 
allel or perpendicular to the surface, represented in red and 
blue, respectively). The full curves show the Bose-enhanced 
contribution EC (JV, d) , while the dashed curves give the 
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Figure 3: Bose-enhanced decay rate ■y e (N) near an interface. The 
BEC contains 10 5 Rb atoms in an oblate wave function at an aver- 
age distance d from a perfectly reflecting surface. The size parameters 
are an = 3.4 (im parallel and an/\/A = 0.01 an perpendicular to the 
surface, corresponding to trapping frequencies of ui\\j2iT = 10 Hz and 



10 u)u and a critical temperature T t 



(2D) 



118 nK. The excited 



wave packet (resonance frequency as in fig. Q is spatially centered in 
the BEC, with size parameters n = 0.07 and rj z = I relative to an. Full 
red curve: Bose-enhanced decay rate 7;? EC (-/V) given by eqn. 13b, for 
an excited atom with its dipole moment oriented parallel to the surface. 
We normalize to the free-space decay rate 7^°' (d — > 00). Full blue curve: 
7? EC (-^0 f° r an excited atom with perpendicular dipole moment. Hor- 
izontal dashed lines: asymptotic values of r y^ EC (N, d — > 00) at large 
separation. Dashed red (blue) curve: single-atom decay rate 7e°' (first 
line of eqn. |3]l) for parallel (perpendicular) dipole orientation; these 
data are multiplied by a factor of 390 and 150, respectively, such that 
their asymptotic values for large distances d coincide with 7 ( ? EC (A r ). 



single- atom part ^^(d), re-scaled such that the asymptotic 
value for large distances d coincides with 7 ( ? EC (iV, d — > 00). 
The numbers given in fig. [3] are the result of a compromise 
between a tight confinement in the vertical (z-) direction and 
a localized wave packet in the excited state. The atomic 
wave packets must be confined below the wavelength in the 
z-direction, otherwise the oscillations in "f e vs. distance are 
averaged out. In this limit, the optimal Bose enhancement 
is found for an excited wave packet that is matched to the 
condensate (?/ z = 1). For the size parameter in the xy- 
plane, we find an optimum at r\ = 0.07. The asymptotic 
values 7 e (iV, d — > oo) are enhanced by factors of 390 and 150 
compared to ji° (d — > 00) for the parallel and perpendicular 
dipole, respectively. The difference between these two num- 
bers and the relative phase shift of the oscillation pattern in 
EC (d) are due to the radiation pattern of the dipole emis- 
sion, combined with the shape of the ground state mode that 
modulates the Bose enhancement in k-space. 

Fig. [3] thus demonstrates a significant amplification of the 
decay rate above the surface, with the oscillation amplitude 
receiving an additional enhancement relative to the asymp- 
totic free-space component. It suggests that even at a dis- 
tance of a few microns (several transition wavelengths) , Bose 
enhancement can bring the tiny interference structure of the 
decay rate into an experimentally detectable regime. 



4 Virtual excited atoms produced by laser absorp- 
tion 

The calculation above assumed the presence of an excited 
atom prepared in a gaussian wave packet, and one may ask 
the question whether this is a realistic description. Indeed, 
the preparation of such a state would typically proceed by 
illuminating the system. We therefore describe in this sec- 
tion a calculation of a typical absorption spectrum. We find 
that the results of the previous section are qualitatively un- 
changed. The method also illustrates the relevance of two- 
and four-point correlation functions of the Bose gas. For the 
sake of simplicity, we restrict this analysis to the ideal Bose 
gas. 

The calculation proceeds by keeping a continuum of modes 
for the excited state field operator ^ e (x) and by identifying 
the absorption spectrum of a weak laser field with a suitable 
T- matrix element (self-energy). We take the laser field to be 
described by a coherent state |/3) in a given plane- wave mode. 

In the leading order of perturbation theory, the absorption 
by the atom cloud of a photon out of the coherent state |/3) 
and re-emitting it into the same state, 



(10) 



results in a (complex) energy shift of the laser plus atom 
system that is described by the T-matrix element 



(N,p\TW\N,f3) = 



e a (k L )ef3(k L ) 



UJ L -IE 



In eqn. ( 11 ) 



(11) 

is the number of photons in the coherent 
state, lol and denote the frequency and wave vector of 
the absorbed photons, the unit vectors e(krj denote axes 
of (linear) photon polarization, and the infinitesimal e \ 
ensures the adiabatic switching-on of the laser field. At this 
order of perturbation theory, the absorption of photons by the 
atom cloud is proportional to Im (TV, f3\T^ \ N, /?) oc 5(ui eg — 

The next order in perturbation theory brings about the 
diagram 

N ^ T . r* ^ . r J -N , (12) 

which gives the following contribution to the T-matrix: 

(N,P\T (4) \N,P) = -l^eMe^)-. ^ 

2 [uj eg - uj l - iey 

x j d 3 Xl J d 3 x 2 (*t( X2 )$ ff ( X2 )^t( Xl )^( Xl )) 

x Al7A/5 G , 75 (x 1 ,x 2 ,w L )e^-( Xl - X2 ) . (13) 

Let us introduce the density correlation function of the Bose 
gas as 

C(x 2)Xl ) = (*t( X2 )^( X2 )^t (xi) ^( Xl)) _„( X2 )„( Xl ) 

(14) 

where n( Xl ) = (ty] ^(x)^ 9 (x)) is the average density. This 



splits eqn. (131 in two parts: (T^) = (T s ( c 4 a } t ) + (T^). The 



former contains only densities and can be identified with the 
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elastic scattering of photons off the inhomogeneous density 
profile of the BEC. This term does not distinguish between 
a Bose gas and a classical system with the same density. Its 
imaginary part provides, by the optical theorem, the total 
scattering cross section of the BEC. The second term (T^ s ), 
on the contrary, is proportional to density fluctuations, and 
these are at the origin of bosonic enhancement [231 HI] ■ We 
therefore identify Im (T^bs) with the change in the atomic 
absorption spectrum (line width). 

Indeed, if we define the resonant part of the polarizability 
a a p(iS) of the atom cloud as 



(15) 



"eg 



we see that the process eqn. ( 12 ) can be re- written as a shift 
of the atomic transition frequency uj eg 



Suj eg with 



2h 



<7i 4 l> = 



e Q (k L )e/3(k L ) 



eg 



(16) 



By identifying eqn. (13) and eqn. (161, we can read off the 
frequency shift Suj eg whose imaginary part yields the atomic 
line width (the inverse lifetime of the virtual state involving 
an excited atom) 

7abs(A0 = -21m Sui eg 



- JjHaUp 



J d 3 xn(x) ImG a p(x, x, lul) 



+ J d 3 Xl y"rf 3 a: 2 C(x 1 ,x 2 )ImG Q(3 (x 1 ,x 2 , WL ) e -^ ( Xl - X2 ) 

(17) 

This function depends weakly on the laser frequency, and we 



by squares). The size parameter 77 was set to the optimal 
value obtained from Fig [I] Although the line width 7 a bs(A0 
is for these parameters around 30% smaller than the opti- 
mized decay rate r y e (N), the strong Bose enhancement is still 
working in a qualitatively similar way for both types of pro- 
cesses. We expect a similar result to hold for an absorption 
experiment near a surface, using for example evanescent fields 
as discussed in Refs. |2l)I fI7\ |2"g] . 

5 Summary 

To summarize, the presence of a trapped BEC can signifi- 
cantly enhance the decay of an excited atom by bosonic stim- 
ulation. The magnitude of the effect depends on the overlap 
between the atomic wave functions and the wavevector of 
the photon involved in the decay. More precisely, our calcu- 
lations based on a quantum field theory of the atom-photon 
interaction illustrate the importance of two- and four-point 
correlation functions of the ground-state field for the Bose 
enhancement. For an excited atom prepared in a Gaussian 
wavepacket, the transition rate to the ground state can be 
increased under optimum conditions by a factor N/10 where 
N is the atom number in the BEC. This effect also amplifies 
the small oscillations of the decay rate near an interface. We 
have provided an alternative calculation based on the absorp- 
tion of a laser beam that qualitatively confirms the simpler 
wave packet picture. The main difference is that absorption 
from the laser field prepares a non-pure excited state which 
matches the one-body density matrix of the Bose gas. 
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evaluate it at ljl 



0J t . 



for simplicity. We shall use below References 



Wick's theorem to evaluate the density correlation function 
[25] . as appropriate for the ideal Bose gas: 



C(x ls x 2 )=|<*t( X3 )^( Xl ))| 



(18) 



Eqs.(17 18) can now be compared to the decay rate 7 e (A0 
of an excited-state wave packet (eqn. ([3])). The first term 
in both expressions is very similar, and we see that the 
laser spectroscopy effectively prepares an excited state den- 
sity profile matched to the condensate density. The second 
terms differ because the laser wave-vector appears explic- 
itly. Also the one-body density matrix for the excited state, 
(e\ , i>l(x2) 1 i' e (xi)\e) in eqn. |2), is replaced by its ground-state 
equivalent (N\^ g (x2)^ g (xi) N) in eqn. (17). The prepared 
wavepacket is hence no longer pure. This makes the tem- 
perature dependence of 7abs(A0 stronger, as can be seen in 
Fig{2] (compare the blue and red curves). The calculation of 
7abs(Af) involves, because of the squared correlation function, 
a double summation over single-particle trap states [15] un- 
der the integral. At zero temperature, the summations can be 
done analytically, at the low temperatures Ta = 0.2, 0.3 T c , 
the double sum could be evaluated numerically (circles in 
fig. [2]), while for Ta/T c > 0.5, the summations can be accu- 
rately replaced by integrations that evaluate faster (denoted 
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